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We study the asymptotic joint distribution of sample space-time covariance estimators of strictly 
stationary random fields. We do this without any marginal or joint distributional assumptions 
other than mild moment and mixing conditions. We consider several situations depending on 
whether the observations are regularly or irregularly spaced and whether one part or the whole 
domain of interest is fixed or increasing. A simulation experiment illustrates the theoretical 
results. 
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1. Introduction 

Let {Z(x) :xeflc M. d }, d>l, be a strictly stationary random field with covariance 
function 

C(k) = cov{Z(x), Z(x + k)}, 

where k denotes an arbitrary lag in M. d . Let A be a set of lags and let m denote the 
cardinality of A. Define G = {C(k) : k <G A} to be the length m vector of covariances at 
lags in A. Let C n (k) denote a sample-based moment or kernel estimator of C(k) based on 
observations in a sequence of increasing index sets D n C D and let G„ = {C„(k) : k £ A} 
denote the estimator of G. We are interested in the asymptotic distribution of G n . This 
distribution is important for several reasons, such as assessing the covariance structure 
of a random field. For example, Li, Genton and Sherman (2007) proposed a testing ap- 
proach for various properties of space-time covariance functions based on the asymptotic 
distribution of covariance estimators in a particular setting. For the particular case of 
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testing directional properties of a random field, Guan, Sherman and Calvin (2004) and 
Lu and Zimmerman (2001) derived the asymptotic distribution of the spatial variogram 
when d = 2. Surprisingly, the distribution of sample covariances has been investigated 
only in particular situations in the literature. Under the assumption of a Poisson process 
for modeling the observations' locations, Masry (1983) proved the asymptotic joint nor- 
mality of sample autocovariances for time series and Karr (1986) generalized this result 
to a random field. Brockwell and Davis (1991), page 229, and Fuller (1996), page 333, 
derived the asymptotic joint normality of sample autocovariances under mild assump- 
tions for stationary time series. For spatially indexed observations, several authors (e.g., 
Bolthausen (1982), Guyon (1995)) have proven the asymptotic normality of the scalar 
sample mean under a variety of mixing and moment conditions. 

Recognizing various situations occurring in practice, we consider several data structures 
and asymptotic regimes, depending on whether the observations are regularly spaced or 
irregularly spaced and whether a subset of dimensions is fixed or the whole domain of 
interest is increasing. We also allow one dimension to denote time. We make no assump- 
tions on the marginal or the joint distribution of observations other than mild moment 
and mixing conditions on the random field. 

To formally state the asymptotic properties of the vector of sample covariances G n , 
we need to quantify the strength of dependence in the random field, taking account of 
different types of spacing of observations. Following Rosenblatt (1956), we define the 
strong mixing coefficients for a random field with regularly spaced observations as 



where \E\ denotes the cardinality of the set E, 3(E) denotes the cr-algebra generated 
by the random variables {Z(x) :x G E} and where d(E\,E2) = infjsupj |xy — X2j| :xi g 
Ei,X2 S E 2 ,j = lj • • • j d}. The supremum in (1.1) is taken over all compact and convex 
subsets Ei C M. d and E 2 C K d such that d(Ei,E 2 ) > r. For a random field with irregularly 
spaced observations, we define the mixing coefficients following Politis, Paparoditis and 
Romano (1998). This definition, denoted (1.1'), is formed by imposing an additional 
condition that E2 is a shift of E\ in (1.1). Note that \E\ denotes the Lebesgue measure 
(volume) of E in (1.1'). 

If the observations are independent, then a&(r) = for all r > 0. We need ctb(r) to 
approach for large r, at some rate depending on b. In order to appropriately define this 
rate, we decompose D n into D n —Txl ni where T <zW,I n C R q and p + q — d. Suppose 
that T is a fixed space, in the sense that finitely many observations arc located within 
this space, and I n is an increasing space. Following Sherman and Carlstein (1994), we 
assume the mixing condition for I n 



a b {r) = sup {\P{A X n A 2 ) - P{A l )P{A 2 )\ : A, e £(£,), \E,\ < 6, 



^1 



(1.1) 



i = l,2,d(E 1 ,E 2 )>r}, 



sup 

b 



ctbjr) 
b 



0(r" s ) 



for some e > q. 



(CI) 



From Proposition 3.1 of Kiinsch (1982), it can be shown that (CI) holds for a class of 
Gibbs fields. Bradley (1993) has also studied and justified the use of mixing condition 
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(CI). We account for the shape of the domain in which we observe data as in Sherman 
(1996). Let A denote the interior of a closed hypersurface contained in a g-dimensional 
hypercube with edge length 1. A n is the inflation of A by a factor n. We define l n = 
1 q D A n if the observations are regularly spaced, and X n = A n otherwise. Note that I n 
satisfies 

\T n \=0{n% \dl n \=0(n«- 1 ), (C2) 

where dl n denotes the boundary of I n . This allows for a wide variety of domains. 

If a random field provides sufficient information so that we can estimate its covariance 
at any arbitrary lag, this random field is said to exhibit continuous lags. We use a 
kernel estimator to estimate covariances in this case. For example, if the observations 
are irregularly spaced in the subspace R 91 of R 9 (<ji < q), covariances for any lag value 
in R 91 can be estimated by a kernel estimator; see Section 3. However, if covariances at 
only discrete lags are estimable, this random field is said to exhibit discrete lags and we 
use a moment estimator to estimate the covariances in this situation. 

Throughout the paper, we assume that the mean of Z is known and equal to 0. If 
we remove this assumption, let C*(k) and G* denote the mean corrected estimators 
of C(k) and G, respectively. In Lemma A. 6, we show that G* and G n have the same 
asymptotic properties. Let C„(k) denote a moment estimator or a kernel estimator of 
the covariance C(k) under the zero-mean assumption and let A„ be the bandwidth of 
the kernel estimator defined over an irregularly spaced subspace R 91 of R 9 . We assume 
the following moment condition for the moment estimator: 

su P E{| v /l2^|{(7 n (k)-C(k)}| 2+ ' 5 }<C5 for some 5 > 0, C s < oo. (C3) 

n 

The moment condition (C3') for the kernel estimator can be obtained from (C3) by 
simply replacing \T n \ with |I„|A 91 and C(k) with E{C„(k)}. The moment condition is 
only slightly stronger than the existence of the (standardized) asymptotic variance of 
C n (k). 

In this article, we derive the asymptotic joint distribution of sample covariances for 
space-time random fields {Z(s, t) : s G R 2 , t £ R}. However, the results for this R 2 x R 
space can be easily extended to the R d (d > 3) space. Note that we consider an increasing- 
domain asymptotic framework. For a recent comparison of infill and increasing-domain 
asymptotics and the consequences for maximum likelihood estimators of covariance pa- 
rameters, see Zhang and Zimmerman (2005). 

The rest of the paper is organized as follows. Section 2 demonstrates the asymptotic 
joint normality of sample space-time covariances for discrete lags, while Section 3 ad- 
dresses the distributional behavior in the case of continuous lags. These two sections con- 
sider space-time data structures which are common in applications. Section 4 presents a 
simulation experiment. Appendix A states and proves some useful lemmas. Appendix B 
contains proofs of all theorems. 
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2. Discrete lags 

2.1. Regularly spaced observations with an increasing 
spatio-temporal domain 

Consider a strictly stationary mean zero spatio-temporal random field {Z(s, t) : s € R 2 , t 6 
R}. Let D n = S n x T n be a finite set of lattice points in Z 2 x Z at which observations are 
taken. We allow the lattice in S n and T n to be defined in different metrics. Let h denote 
a lag in space and u denote a lag in time. The classical estimator of the covariance, that 
is, the sample covariance, is given by 



where the sum is over D n (h,u) = {(s,t) : (s,t) £ D n , (s + h.t + u) g D n } and \D n (h, u)\ 
denotes the number of distinct elements in D n (h,u). 

Wc define the strong mixing coefficients of this random field as in (1.1) and assume 
the mixing condition, boundary condition and moment condition as in (CI), (C2) and 
(C3) setting p = 0, q = 3 and T n = D n . 

Theorem 1. Let {Z(s,():s£l 2 ,t£l} be a strictly stationary mean zero spatio- 
temporal random field observed at lattice points in D n C Z satisfying condition (C2). 
Assume that 



for all finite hi, h2,ui and u^. Then, £ = linin^oo \D n \ cov(G„, G„) exists, the (i,j)th 
element of which is 



If we further assume that E is positive definite and that conditions (CI) and (C3) hold, 
then y/\D n \(G n - G) N m (0, S) as n -> oo. 

2.2. Observations with a fixed spatial domain and an increasing 
temporal domain 

In many situations, the observations are taken from a fixed space S CM 2 at regularly 
spaced times T n . Let 5(h) = {s : s € S, s + h E S} and |5(h)| be the number of elements 
in 5*(h). In this situation, we have p = 2, q = 1, T = S and X n = T n in the notation of 
the Introduction. 





(2.1) 



cov { z (°' 0)Z{hi, Ui ), Z(s, t)Z(s + hj, t + uj)}. 



sez 2 tez 
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In this particular case, we define the mixing coefficient (e.g., Ibragimov and Linnik 
(1971), page 306) only in the time direction since the domain is increasing only in this 
direction: 

a(u) = sup{|P(A n B) - P(A)P{B)\,A G^Be 

A.B 

where B .^ is the a-algebra generated by the past-time process until t = and is 
the cr-algebra generated by the future-time process from t = u. We assume the mixing 
coefficient a(u) satisfies the strong mixing condition: 

a{u) = 0(u~ E ) for some e > 0, (2.2) 

and we also observe that 

T n = {l,...,n}, |9T„| = 2 = 0(1). (2.3) 

Condition (2.3) is a special case of (C2), while condition (2.2) holds for a variety of 
temporal processes. For example, (2.2) holds for AR(1) processes with normal, double 
exponential or Cauchy errors by the results of Gastwirth and Rubin (1975). We further 
assume the moment condition as in (C3) and define the estimator of C as 

^ n—u 

dn ( h ' u ) = io^ikt i E E z ^ *)^ s + ht + 

l&W||i„| 5(h) t=1 



Theorem 2. Let {Z(s,t),s £ M. 2 ,t £ K} be a strictly stationary mean zero spatio- 
temporal random field observed in D n = S x T n , where ScK 2 and T n satisfies condition 
(2.3). Assume that 

I cov{Z(0, 0)Z(h u ui), Z(s, t)Z(s + h 2 ,t + u 2 )}\ < oo (2.4) 

tez 

for all hi <= S, h 2 € S, s £ iS(ha) and all finite u\,u 2 . Then, S = linin^oo \T n \ cov(G„, G„) 
exists, the (i,j)th element of which is 

lev, -} lQf u \\ E E 'V.cov{Z(s 1 ,0)Z(s 1 + h i ,Ui),Z(s2,t)Z(s2 + h. j ,t + u j )}. 

If we further assume that £ is positive definite and that conditions (2.2) and (C3) hold, 
then \/\T n | (G n — G) — >N m (0, S) asn^oo. 



In Theorem 2, we allow the observations to be either regularly spaced or irregularly 
spaced in S. However, even for irregularly spaced observations, we consider only the 
covariances of observed spatial lags due to the limited number of observations in S. Note 
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that in this section, we require the observations to be taken at the same spatial locations 
over time, which is very common for monitoring stations. For example, the Irish wind 
data of Haslett and Raftery (1989) recently analyzed by Gneiting (2002), de Luna and 
Genton (2005) and Stein (2005) consists of time series of daily average wind speed at 
eleven meteorological stations in Ireland. In particular, Li, Genton and Sherman (2007) 
developed space-time symmetry and separability tests based on the asymptotic normality 
of sample covariance estimators under this specific data structure. 



3. Continuous lags 

3.1. Spatially irregularly spaced observations with an increasing 
spatio-temporal domain 

If the observations are spatially irregularly spaced in an increasing domain, we can then 
estimate the covariance for any spatial lag by employing kernel smoothing. Again, con- 
sider a strictly stationary mean zero random field {Z(s,t) :s <G R 2 ,£ € K}. Let D n C R 3 
denote the domain of interest in which observations are taken. We decompose D n into 
D„ = S„ x T n , where S n is an increasing index set, S n C R 2 and T n = {1, . . . , n}. We view 
the spatial locations at which Z is observed in S n as random in number and location; 
specifically, they are generated from a homogeneous 2-dimcnsional Poisson process with 
intensity parameter v. 

Let N denote the random point process and N(B) denote the random number of 
points of N contained in B, where B is any given Borel set. We further assume N to be 
independent of Z . We construct an estimator of covariance based on kernel smoothing. 

Let \S n \ denote the Lebesgue measure (not the cardinality) of S n and let w(-) be a 
bounded, symmetric density function on R 2 . Here and henceforth, we use ds to denote an 
infinitcsimally small disc centered at s. Define A r ^ 2 - ) (dsi,ds2) = -/V(dsi)7V(ds2)/(si ^ S2), 
where I(si ^ S2) = 1 if Si ^ S2 and otherwise. The kernel covariance estimator over D n 
is given by 

^ n—u ~ ~ 

C n (h,u) = -^p || g 1 J2 J s J s ^r l (h-s 1 +s 2 )Z(s 1 ,t)Z(s 2 ,i + u)^ (2) (ds 1 ,ds2), 

where w n (x) = -p-w(-^-) and A„ is a sequence of positive constants satisfying condition 
(3.1) below. Here, v is assumed to be known, otherwise it can be consistently estimated 

, - N(S n ) 

We define the mixing coefficients as in (1.1') and assume the mixing and moment 
conditions as in (CI) and (C3'), setting p = 0, q = 3, q\ = 2 and I n = D n . In addition to 
the boundary condition in (C2) , we need to consider the choice of bandwidth. Specifically, 
we assume that 

A„^0, X 2 n \S n \^oo. (3.1) 
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Define the following fourth order cumulant function (e.g., Karr (1986)): 

Q(x 1; x 2 ,x 3 ) = E{Z(0)Z( Xl )Z(x 2 )Z(x 3 )} - C( Xl )C(x 3 - x 2 ) 

(3.2) 

- C(x 2 )C(x 3 - xi) - C(x 3 )C(x 2 - Xl ), 

where the x's denote generic locations and x = (s, t) in the spatio-temporal random field. 

The following theorem states that C n (h,u) is a consistent estimator of C(h,u) and 
G„ is asymptotically jointly normal under some conditions. 

Theorem 3. Let {Z(s,i):seR 2 , f £ K} be a strictly stationary mean zero spatio- 
temporal random field observed in D n = S n x T n , where S n C K 2 and T n = {1, . . . , n}. 
D n satisfies condition (C2). Assume that the locations s are generated by a homogeneous 
Poisson process in M. 2 and that 

sup / C(h,it)dh < oo, 

sup '^E{Z(O,0)Z(h,u)Z{O,t)Z{h,t + u / )} < oo, 
h,u,u' tgz 

sup sup / |Q{(si,ti), (s 2 ,t 2 ), (s 2 + s 3 ,i 2 +t 3 )}|ds 2 < oo. 

SI, S 3 tl,t 2 ,t 3 J S2 gR2 



s 2 Gl 

v2 



Then, E{C„(h, u)} — > C(h, u) and T, = lim n ^ OD \T n \\S n \Xn cov(G n ,G n ) exists, the (i,j)th 
element of which is 

4 / w 2 ( X )d X yE[Z(0,0)Z(h t ,u l ){Z(0,t)Z(h l ,t + u j )L(h i =h J ) 



+ Z(h i ,t)Z(0,t + u j )I(hi = -h j )}], 

where L(u = v) = 1 if u = v and otherwise. If we further assume that E is positive 
definite and the conditions (CI), (C3') and (3.1) hold, then 



^J\T n \\S n \\ n {G n - E(G„)} N m (0, E) 



as n — ■> oo . 



3.2. Irregularly spaced observations with an increasing 
spatio-temporal domain 

In Section 3.1, we discussed the properties of the kernel covariance estimator if the 
observations are irregularly spaced in S n and regularly spaced in T n . Another case occurs 
when the observations are irregularly spaced in the whole space S n x T n . Karr (1986) 
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showed mean square consistency of the estimator C of the covariance function C of the 
random field Z defined on R d (d > 2). It is an extension of the results in Masry (1983) 
who investigated the covariance estimator for time scries, that is, for random fields in 
R 1 . 

In Karr's theorem, the random process, {Z(pt) :xG R d }, is assumed to be a stationary 
Poisson process with intensity v which is independent of the values of Z(x). We specialize 
his results to our space-time random field by considering {Z(x) :xgl 3 }. Let x denote 
a location in R 3 and dx denote an infinitesimally small sphere centered at x. Define w(-) 
and ./V( 2 - ) (dxi,dx2) in the same way as in Section 3.1 with the understanding that in this 
section, the support of w(-) and x are defined in R 3 rather than R 2 . The kernel estimator 
over D n is defined as 

C^) = ^Wli i «>n(k-x 1 +x 2 )Z(x 1 )Z(x 2 )iV< 2 >(dx 1 ,dx 2 ). 
v*\D n \ J Dn J Dn 

We define the mixing coefficients as in (1.1') and assume the mixing condition, the 
boundary condition and the moment conditions (CI), (C2) and (C3'), setting p = 0,q = 
3, qi = 3 and X n = D n . We adopt the definition of Q from (3.2) in this section and assume 
that A„ satisfies the condition 

A.,,^0, X 3 n \D n \^oo. (3.3) 



Theorem 4. Let {Z(x):x£l 3 } be a strictly stationary mean zero spatio-temporal 
random field which is observed in D n C R 3 satisfying condition (C2). Assume that 
the locations x are generated by a homogeneous Poisson process in R 3 . Assume that 
J R3 C(k) dk < oo and that Q exists and satisfies 

sup / |Q(x + xi,x,x 2 )| dx < oo. 

xi,x 2 JK 3 

Then, E{C(k)} — > C(k) and S = linin^oo A 3 \D n \ cov(G n , G„) exists, the (i,j)th element 
of which is 

-1 / ^ 2 (x)dx[E{Z 2 (0)Z 2 (k i )}/(k l =±k,)], 

where /(k, = ±kj) = 1 if k^ = ±kj and otherwise. If we further assume that 
E{Z 2 (0)Z 2 (k)} > for alike A and that conditions (CI), (C3') and (3.3) hold, then 

VA 3 |£> n |{G„ - E(G„)} N m (0, E) 

as n—too. 

Unlike Theorem 3, Theorem 4 combines s and t as a single location x. As a result, we 
get a more concise form for S than in Theorem 3. However, we can expect there to be 
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two separate terms for /(hi) = I(hj) and 7(hj) = I(—h.j) in the formulation of £ as in 
Theorem 3 if we consider s and t separately. 

4. Simulation 

In order to illustrate the approach to joint normality of G„ with the derived asymptotic 
covariance matrix, we perform the following simulation experiment. For simplicity, we 
consider only the situation given in Section 2.2, noting that it corresponds to many 
practical applications such as the Irish wind data mentioned at the end of Section 2.2. We 
simulate a random field via a vector autoregressive (VAR) model with spatial structure 
(see, e.g., de Luna and Genton (2005)). Specifically, we assume that the fixed space S is 
a 3 x 3 grid with grid interval 1 and consider the VAR(l) model 

Z t =RZ t _ l +e U (4.1) 

where Z t = (Z(s\,t), . . . , Z(sg, t)) T and is a Gaussian multivariate white noise process 
with a spatial stationary and isotropic exponential correlation function given by (£ e )jj = 
exp(— H s '~ s ?ll ), i, j = 1, ... ,9, where <j> denotes the range parameter (we initially set <f> = 
1). R is a 9 x 9 matrix of coefficients which determines the dependency between Z t and 
Z t _ r as cov(Z t , Z t -r) = i? r I\(0), where vec-jT 2 (0)} = (I 81 - R® i?)" 1 vec(£ e ) and vec(-) 
denotes the operator vectorizing a matrix. Note that the coefficients in R are determined 
only by s, not by t. We set these coefficients as follows: for each (s,,t), i = 1, . . . , 9, the 
coefficient is 0.2 for (sj,i — 1), whereas it is 0.1 for {(Sj, t — 1) : \\sj — Sj|| = 1}, 1 < j < 9, 
and for the remaining (s,t — l)'s. 

We choose two space-time lags, ki = (||h|| = 1, u = 0) and k2 = (||h|| = 1, u = 1), and set 
A = {ki,k2}. Here, |jh|| denotes the Euclidean norm of h. To explore the empirical joint 
distribution of G„ as T n = {1, . . . , n} increases, we first set n = 3, simulate an5x T n ran- 
dom field using (4.1) and compute C n (ki) and C , „(k 2 ) over this random field. We repeat 
this 5000 times to obtain a sample of (7 ra (ki) and C„(k2) values and then compute their 
empirical covariance matrix. We then set n = 10, 20, 50, 70, 100, 150, 200, 500, 1000, 5000, 
respectively, and repeat the same procedure as with n = 3. We assess the joint normal- 
ity of C n (ki) and C„(k2) for each n using the multivariate measures of skewness and 
kurtosis introduced by Mardia (1970). 

Finally, we evaluate the covariance matrix of G„. Denote by Cij(u) the sample estima- 
tor of Cij(u) = cov{Z(sj,t), Z(si,t + u)}. Priestley ((1981), page 693) presents a formula 
to compute cov{Ci 1 j 1 (ui), Ci 2 ,j 2 (u2)} for large n and Gaussian processes. The result of 
this is essentially an immediate form of Theorem 2 once we assume that Z is a Gaussian 
random field and we are solely interested in Cij(u) rather than C„(k). However, to ob- 
tain the theoretical values for the general C n (k) in Theorem 2 using Priestley's method, 
we need to find a link between the form in Theorem 2 and Cij(u). 
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n 


bi,2 


62,2 


|T n |cov{C n (ki),C n (k 2 )} 


|T„|var{C„(ki)} 


|Tn|vax{C„(k 2 )} 


3 


5.797 


17.565 


0.382 


0.491 


0.559 


10 


2.730 


12.929 


0.510 


0.651 


0.584 


20 


1.291 


10.279 


0.483 


0.631 


0.542 


50 


0.701 


9.342 


0.497 


0.643 


0.553 


70 


0.435 


8.917 


0.498 


0.637 


0.552 


100 


0.270 


8.583 


0.492 


0.647 


0.533 


150 


0.185 


8.447 


0.505 


0.663 


0.546 


200 


0.179 


8.352 


0.494 


0.652 


0.533 


500 


0.056 


8.026 


0.496 


0.655 


0.538 


1000 


0.043 


8.084 


0.511 


0.665 


0.554 


5000 


0.010 


8.053 


0.501 


0.653 


0.539 


oo 





8 


0.497* 


0.653* 


0.539* 



61,2 and f>2,2 denote multivariate measures of skewness and kurtosis, respectively. 
*are obtained using Lemma 1. 

Lemma 1. If Tit = (Z(si, t), . . . , Z(s n , t)) T is a Gaussian process with mean 0, Cij{u) = 
E{Z(sj, t)Z(si, t + u)} and Cij(u) = yjr-^^2 t Z(sj,t)Z(si,t + u), then 

cov{C n (h i ,u i ),C n (h j ,u j )} = —^^—-^ co^{^k,k'(ui),Ci i i l (u j )}, 
where sy = Sfe + hj, Si> = si + hj, and for large n, 

cav{d iuh (s),C i2th (u)} - — ^{C i2;il (r)C l2 . M (r + u~s) + C i2tjl (r + u)C 32 . n (r - s)}. 

We calculate the theoretical values in Theorem 2 using Lemma 1 and compare them 
with the simulation output. The simulation output, together with the theoretical val- 
ues, are summarized in Table 1. We see from this table that the multivariate measures 
of skewness and kurtosis approach and 8 respectively as n increases. This agrees 
with the skewness and kurtosis of a bivariate normal distribution. Additionally, all 
|T„|cov{C n (ki),C„(k 2 )}, |T„| var{C„(ki)} and \T n \ var{C„(k 2 )} approach their corre- 
sponding theoretical values as n increases. This verifies our asymptotic covariance matrix 
of G„. As seen in Table 1, cov(G„, G„) stabilizes at about n = 100. 

We have observed that the simulation results do not change appreciably when the 
number of simulations exceeds 5000, so we show only the results based on 5000 simula- 
tions. Setting the range parameter </> in the spatial correlation function to be 1.5 increases 
the variances and covariances, but the trends of multivariate measures of skewness and 
kurtosis and the convergence of cov(G„, G n ) remain the same as when = 1. When S is 
a 4 x 4 grid, the variances and covariances become smaller than with a 3 x 3 grid as larger 
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Figure 1. Partition of the random field for Lemma A. 2. 

grid sizes provide more data to estimate G„. However, all of our simulation experiments 
show a more rapid approach of G„ to joint normality as n increases. 

Appendix A: Lemmas 

Lemma A.l. Consider two closed and connected sets U , V in M. d such that \U\ = \V\ < b 
and d(U, V) > r. Let X and Y be measurable random variables with respect to 3(U) and 
$(V) such that \X\ < Ci and \Y\ < C 2 . Then, cov(X,Y) < ^CiC 2 a b (r). 

Proof. The proof is completely analogous to that of Theorem 17.2.1 in Ibragimov and 
Linnik (1971), page 306. Note that if the variables X, Y are complex, then separating 
the real and imaginary parts, we again arrive at the same expression, with 4 replaced by 
16. □ 
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Lemma A. 2. Let {Z(x),x £ M. d } be a strictly stationary mean zero random field which is 



Then, a 2 =lim n _ f00 \D n \Yai(Tr n ) exists. If we further assume conditions (CI), (C2) and 
(C3), setting p = 0, q = d, X n = D n , C„(k) = 7?„ and C(k) = 7r, then y/\D n \(n n — n) — — > 



Proof. Let A„ = \/\D n \(w n - tt). Then, a 2 = lim, woo var(A„) = J2xez d cov [/{^(°)}> 
/{Z(x)}] exists by assumption. We apply a blocking technique (e.g., Ibragimov and 
Linnik (1971)) in conjunction with the mixing condition to prove the normality of 7? n . 

Let l(n) = n a and let m(n) = n a — n v for some 2d/(d + e) < r\ < a < 1. Divide the 
original field D n into non-overlapping subhypercubes, DJ,-, =l d (n), i = 1, . . . , fc n ; within 
each subhypercube, further obtain D l m ^ which shares the same center as DJ, n y Thus, 
d{D l , n y D l m i n \) > n v for i^i' . Figure 1 shows an example when d = 3. Let TT l m r n \ denote 
the estimator obtained from D % m i n \ ■ Let a n = J2i=i a nl V^n, a ' n = Yli=i ( a nY / Vkn, where 
a' n = \fm d {n)\Tr' l m< s — tt} and (a % n )' have the same marginal distributions as a l n , but 
are independent. Let 4>' n {x) and 4> n (x) be the characteristic functions of a' n and a n , 
respectively. The proof consists of the following three steps: 



(51) A n -a n AO; 

(52) <(/ n (x)-M*)^0; 

(53) <4N(0,a 2 ). 

Proof of (SI). Since E(A n — a n ) = 0, it suffices to show that var(A n — a n ) — > 0. Let 

D m(n) denQte the 

union of all D*,^. Observe that 




]T |cov[/{Z(0)},/{Z(x)}]|<oo. 



N(Q,a 2 ). 



1 



A: 




n 





|£)m(n)| 
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Thus. var(a„) — > a 1 since rp^tyi — ► 1 (demonstrated in Lemma A. 3) and 



/In \\r)m(n)\ 

^■'"^ Ifltw £ £ cov[/{Z(x 1 )},/{Z(x 2 )}]. 

Again by . I — > 1, we get cov(yl n , a n ) — > <r 2 and var(A„ — a„) — > 0. 

Proof of (S2). We employ telescope arguments here. Let i denote the imaginary unit. 
We define Ui = exp(tx^=), Xj = Yli=i Ui and Yj — Uj+\- Analogously to Sherman and 
Carlstein (1994), Lemma 2, we have 

covp^-, Y,) < lGjm d (n)n- er i = \Qj(n a - n r ') d n _£ '' 

by (CI). Thus, \cf>' n (x) - <j) n (x)\ < 16jn da - En = 0(n 2d - da - E7 >) . The last equality 

in the previous expression follows from 0(k n ) = 0(~^) = 0(n d ( 1_Q )), and X^" 1 ^3 = 
0(k 2 l ), Since 2d/(d + e) < 77 < a < 1, we have 2d - da - £77 < 2d - drj - £77 < 0. Then, 

Proof of (S3). Observe that E(|(a^)'| 2+(5 ) < C& for some constant C$. Since [a l n )' are 
i.i.d., var{^ i ^ 1 (a^ l )'} = k n var{(a^)'}. Defining a 2 = var{(a^)'}, we have <r 2 — > a 2 from 
the proof of (SI). Thus, 

„ m £ gflMffla < lim c *. p. 

i-.oo^-' /; r „fr„ „ , c ~~ ra->oo k„a 2 W+ s )/ 2 



n — kx> 

i=l 



A rv- fe " / » \/T12+5 ~ (fcn^) (2+<5)/2 



Thus, applying Lyapounov's theorem, we have -7= X)i=i( a n)' — >N(0,er 2 ). □ 

Lemma A. 3. Assuming that D n satisfies (C2) setting T n = D n and q = d, we have 
\D m ( n '\/\D n \ — > 1 as n—> 00, where D m (") is defined as in the proof of Lemma A. 2. 

Proof. Introduce the following notations: D'W denotes the union of -D|(„) defined in 
the proof of Lemma A. 2; D n /u n ) denotes the field in D n but not in D 1 ^; Dii n )/ m (n) 
denotes the field in D 1 ^ but not in D" 1 ^; k'[n) denotes the minimal number of extra 
l d (n) subhypercubes needed to cover the whole of D n ; D n i denotes the union of all extra 
l d (n) subhypercubes needed to cover the whole of D n . Since 

\D l W\ = K\D\ (n) \ = k n l d (n) = k n n da , 
\D m ^\ = k n \D l m[n) \ = k n m d (n) = k n (n a - n r ') d = k n n da + o{k n n da ), 
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noting that V < a, we have |A(»)/m(»)l = \D l{n) \ - |^ m(n) | is o(\D L ^\). 

We now show that |-D n / i(n) | is o(n d ). Note that D n/Hn) C.D n ', that is. |D n /j (n )| < \D n ,\. 
Thus. \D n t\ being o(n d ) is sufficient for |-D„//(„)| to be o(n d ). 

To show this, we split the boundary of D n into hypercubes of volume Z d_1 (n). Here, we 
use L ni , to denote the ith hypercubc and k',' n to denote the number of all the hypercubes 
available, k^ is 0( l J l _ i(„) ) due to (C2). For each L n ^, we may construct an l d {n) hyper- 
cube (denoted by HCi^ i) which fully contains L n j; in addition to that, we construct 
a {3l(n)} d hypercube (denoted by HC^ n ) ti ) which has the same center as HC^ n ^i. 

Any l d (n) hypercube that intersects L n ^ (and thus intersects HCi( n )^) will be fully 
contained in HCzi( n ),i- Since every subhypercube in D n > intersects the boundary of D n 
and thus at least one of these L nj! 's, it will be fully contained in one of the HC^ n )/s. 
These l d {n) hypercubes in D„/ do not intersect each other, except at the boundary. 
The maximum number of such hypercubes that intersect L n> i cannot be larger than 3 d . 
Since the size of fc" is 0( ;d "'' 1 ,\ ), we conclude that k' n is no larger than 0( l J l _ 1 ,, ) due 



to k' n < 3 d C Thus, \D n/l[n) \ is no larger than 0{^^l d {n)) = 0(n d -H(n)) = o(n d ). 
Then, 



Lemma A. 4. Let {Z(x),x£ M. d } be a strictly stationary mean zero random field which 
is observed in D n = S n x T n , where S n C W , T n C IP and p + q = d. Assume that the 
locations of observations in T n are regularly spaced and that the locations in S n are 
generated by a homogeneous Poisson process. Let \S n \ denote the volume of S n and \T n \ 
the cardinality ofT n . Let 

<?n( k ) = ,,2i T lie? I X! / / w n (kp-Xip + X2p)Z(xip,x g )Z(x 2 p,x 9 +k g ) 




Thus, \D m W\/\D n \ -> 1 as n^oo. 



□ 



x,er, 



where k = (k, 




) T and x = (x, 



x A (2) (dxi p ,dx 2p ), 
:J,x^) T . Assume that 



sup / C(kp, k 9 ) dkj, < co, 

k, JRp 



sup E{Z(0,0)Z(x p , X(Z )Z(0,x l9 )Z(: 

Xf " x " X2 'xi,ez« 



x p ,xi 9 +x 29 )} < 00, 



sup sup / 

xi p ,x 2p x,,xi,,x 2 , Jx p £Ri> 



|Q{(xi p ,xi g ), (x p ,x g ), (xp + x 2p ), (x g + x 2g )}| dx 2p < oo, 
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where w n (x) = jp-w(-^-) and X n is a sequence of positive constants satisfying A„ — > and 
A^|S„| — > oo. Q is defined as in (3.2). Then, E{C„(k)} — ► C(k) and £ = lim rwoo \T n \ x 
|£f n |Af| cov(G„, G„) exists, the (i,j)th element of which is 

^ / w2 (y) d y Hz(o,o)z(k ip ,k iq ){z{o,* q )z{k ip , Xq + k jq )i{k ip = k w ) 

+ Z(k ip ,x 9 )Z(0,x 9 + k jq )I(k ip = -kj p )}], 

where I(ki p = ±kj p ) = 1 if ki p = ±kj P and otherwise. If we further assume that £ 
is positive definite and that conditions (CI), (C2) and (C3') hold, setting p^O, q = d, 
qi =p, I n = D n and replacing C(k) with E{C„(k)} ; then 



T n ||S n |AE{G„ - E(G„)} N rn (0, £). 



Proof. 

E{5„(k)} 



u 2\t\\ s | Y2 E|yy'w„(kp+xi p -X2p)2'(xi p ,x g )Z(x2p,x g + k g )iV (2) (dxip,dx2p)| 
i] E ill// _L W pp ± Xi p X2 p J C(x 2p - xi p , k g ) dx lp dx 2p 



\T n \ ~ \S n \ J J Xn \ X 

1 V- f i \n(i \ i J<S'»n(5'„-k p + A»y)| 
= Ttm 2^ / w(y)C(k p -A„y,k 9 ) — dy, 

where 5„ — SVi denotes all pairwise differences between the two sets. Since 
/ u;(y)C(kp-A n y,k g ) —f dy 

— ►C(k p> k,) / w{y)&y = C{k p ,k q ), 

we have E{C*„(k)} — ► C(k). 

The derivation for the variance is analogous to Karr (1986). Specifically, consider two 
spatial locations k = (k p , k^) T and k' = (k' p , k' q ) T : 

E{5„(k)C„(k')} 

= y4|r n |2 z2 zZ js^2 JJJJ Wn ^ kp + Xlp ~ **p) Wn ( k P + x 'ip ~ x 2p) 

x,GT„ x' GT„ <-, 
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x E{Z(x lp ,Xg)Z(x 2p ,Xg + kg)Z(x' lp ,x g ) (A.l) 

x E{7V< 2 ) (dx lp , dx 2p )7V( 2 ) (dx' lp , dx' 2p )}. 

Observe that 

E{ Z(x lp , x 9 ) Z (x 2p , x g + k q ) Z (x' lp , x' q ) Z (x' 2p , x q + k' q ) } 

= Q{(x 2p -xi p ,k g ), (x' lp - xj p ,x ? - x 9 ), (x' 2p - x lp ,x' q -x q +k' q )} 

+ C(x 2p - x lp , k g )C(x 2p - x' lp , k' q ) (A.2) 
+ C(x' lp — xip, x' q — x g )C(x 2p — x 2p , x ? — Xg + k g — kg) 
+ C( X 2 P - x ip, x g - Xg + k g )C(x 2 p - x' lp ,Xg - x q + kg). 

In addition, 

E{7V< 2 ) (dxip, dx 2p )iV( 2 ) (dx' lp , dx 2p )} 

= is 4 dx ip dx 2p dx' lp dx' 2p 

+ dxip dx 2p er Xlp (dx' lp ) dx' 2p + dx ip dx 2p dx' lp e Xlp (dx' 2p ) (A. 3) 

+ v 3, dx ip dx 2p e X2p (dx' lp ) dx' 2p + ^ 3 dx ip dx 2p dx' lp e X2p (dx' 2p ) 

+ v 2 dx ip dx 2p e Xlp (dx 2p )e Xlp (dx 2p ) + v 2 dx ip dx 2p e X2p (dx' lp )e X2p (dx' lp ) , 

where £ x (") denotes point measure. 

Substituting these two expansions into the main formula produces a lengthy expression 
that we do not reproduce in its entirety. We only show some representative terms, rather 
than the total 28 terms. For simplicity, wc first look at only the integral part conditional 
on Xg and x g . Expanding (A.l), first by (A.2) and then by (A. 3), the first term is 



1 



^\S n \ 



w n (k p + xip - x 2p )w„(k p + x'lp - x' 2p ) 

s„ 

X Q{(X 2p — Xl p , kg), (X lp — Xl p , Xg — Xg), 

( X 2 P ~ x ip , x'g - Xg + k' q )}v A dx ip dx 2p dx' lp dx' 2p 



< // / w n (k p - vi)w„(k' + v 2 - v 3 ) 



S n -S„ 

x |Q{(vi, kg), (v 2 , Xg - Xg), (v 3 , Xg - Xg + k'g)}| T^-f dv i dv 2 dv 3 



1 

1^1 



< Ci J J w n (k p - Vi)w n (kp - v 4 ) |^-r dvi dv 4 
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The second term is 
1 



v 4 \S n \ 



w n (k p + xi p - x 2p )w rl (k / p + x' lp - x' 2p ) 



x C(x 2p - x lp , k,)C(x2 p - x' lp , k' g )j/ 4 dx lp dx 2p dx' lp dx' 5 



2;.» 



/ w n (k p - vi)C(vi,k 9 ) — ^dvi 

JS n -S„ Pnl 

X Js s Wn ^ k p~ V2 ^ C ^ 2 ' k ^ Js~\ 



->C(k) / M (y)dyxC(k') / ™(y)dy 

= C(k)| x ,xC(k')| x; , 
where C(k)| x denotes C(k) conditional on x g . Note that 

W E E C(k)kxC(k')| xi -,C(k)C(k'). 

1 n| x,et„x;gt„ 

Most of the terms arc smaller than or the same as 0( j^-f)- Only 8 among the 28 terms 
are larger than 0(j^-| ) and thus considered as dominant contributions. Combining those 
8 terms, we obtain 

|S„|A£cov{C(k)| x<! ,C(k')|x'J 

— >-W W 2 (y)dyE[Z(0,0)Z(k p ,k ? ) 

x {z(o, x ; - x 9 )z(k p , 4 - x 9 + k;)/(k p = k p ) 
+ z(kp, x ; - Xg )z(o, x ; - x 9 + k;)/(k p = -k;)}]. 

Let 

X « : =K^ E E E[Z(0,0)Z(k p ,k g ) 
|Jn| x,eT„x;eT„ 

X {Z(0,X^ -Xq)Z(k p ,Xq -x q + k' q )I(k p = k' p ) 

+ Z(k pi x.' q - x q )Z(0, x q -x q + k' q )I(k p = -k' p )}} 



E 



Z(0,0)Z(k p ,k,) {^(0,v)Z(k p ,v + k' g )/(k p = k;) 

T„-T„ 
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+ z(k p ,v)z(o,v + k;)/(k p = -k;)} r " n | ^j 2 v) 

Applying Kronecker's lemma, we have 

\T n \X n ^ E[^(0)Z(k){Z(0,x g )Z(k p ,x g + k;)/(k p =k;) 

+ Z(k p ,x,)Z(0,x 9 + k£)J(k„ = -k' p )}]. 

Thus, 

|T„||S n |A£cov{C n (k),C n (k')} 

— J2 E[Z(0)Z(k){Z(0,x 9 )Z(k p ,x 9+ k;)/(k p =k p ) 

+ Z(k p ,x 9 )Z(0,x, +k',)/(kp = -k' p )}] 

x ~2 I u,2 (y) d y- 

v JRp 

There are two different terms for I(k p = k') and I(k p = — ki) because C(k p ,k 9 ) does 
not equal C(— k p ,k 9 ) in general. 

The proof of normality again uses three steps, as in Lemma A. 2. The proof of (SI) is 
analogous to that in Lemma A. 2, but uses a kernel estimator and an additional require- 
ment for a, Xf l n ap — »oo. Wc follow Politis, Paparoditis and Romano (1998) to prove 
(S2). Define Xi and Yi as in Lemma A. 2 and 

E Ar (A J )=E(A l |AT), E N (Y i )=E(Y i \N), cov N (X i ,Y i )=cav(X i ,Y i \N). 

Then, 

cov(A 4 , Yi) = E{cov JV (A J , Yi)} + cov{E w (A 4 ), E N (Y t )}. 

Since N is a homogeneous Poisson process and Xi, Yi are random variables defined 
on two disjoint random fields, we have cov {F> N (Xi),E N (Yi)} = 0. For given N, Xi\N 
is measurable with respect to 3tUj=i ^m(n)) and Yi\N is measurable with respect to 
$( D mL))- We then have WrX x ) - ^ x )\ - cov(X u Yi\N). The rest of the proofs of (S2) 
and (S3) are completely analogous to Lemma A. 2. The proof of the joint normality follows 
from the Cramer-Wold device. □ 

Lemma A. 5. Let {Z(x), x £ R d } be a strictly stationary random field with mean \i which 
is observed in D n C M. d satisfying conditions (CI) and (C2), setting q = d and T n = D n . 
LetZ n ^ J ^- l Z Dn Z( X ). If 



| cov{Z(0), Z(x)}| < oo for all finite x 

x6R d 
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and 

supE| y/\D^\(Z - fx)\ 2+5 < Os, for some 5 > 0, C S < oo, 



then = lim^^cx) \D n \ var(Z„) exists and \j\D n \(Z n — fi) — ► N(0, o~-j^) as oo. 

Proof. The proof follows directly from Lemma A. 2. □ 

It is straightforward to prove that the results in Lemma A. 5 hold for all of the situations 
we have considered in this article, given the appropriate mixing condition. We have proven 
that G„ = {C„(h, u) : (h, u) £ A} is asymptotically jointly normal, that is, \/|-Dra| x (G ra — 

G) — — > N m (0,S) under appropriate assumptions. In this formulation of G n , we have 
assumed the mean of Z(s,t), fJ-(Z), is known and, without loss of generality, assumed 
it to equal 0. However, in most cases, the mean of the random process is unknown and 
needs to be estimated. 

Lemma A. 6. Let ju„ := Z n = pyj X^xec Z(x) an d ^ = {C*(k) :k G A}, where 

Q(k) = \d„{)l)\ ExeD„(k){ 2 '( x ) - Mn}{^(x + k) - /!„}. The vector o/ sampie comri- 

ances G* i/iera has the same asymptotic properties as G n , that is, ^J\D n \{G* n — G) — — ► 
N m (0,E). 



Proof. Assume that fi^O. We now prove that y|D^{C n (k) — C* (k)} = o p (l) as in 
oo. Following Brockwell and Davis (1991), page 229, we have 

v^{c (1 (k)-a;(k)} = v ^— i— Yl (Zn-v){z(x) + z( x + k)} + (v 2 -zl) 

' n \ " x££>„(k) 

= ^\L\\(Z n -^)—^— {%) + Z(x + k)-/i-I„} 



= \J\D n \(Z n - n)(Z n - fi). 



By Lemma A. 5, we know that y/\D n \(Z n — /i) — > N(0,ct-|-), which implies that 
- M) is O p (l) and Z„ - n 0. 
From these observations, we conclude that V /|A^[(G„ -G*) = o p (l). Thus, VlAJ(G* - 
G)-^N m (0,E). □ 

Proof of Lemma 1. When Zt is Gaussian, the fourth order cumulant function Q = 
(Isserlis (1918)). Hence, 

cov{C iltjl (s),C i2tj2 (u)} 

= E{C iujl (s)C i2 , j2 («)} - E{C ilJ - 1 (s)}E{d i2 , j2 («)} 
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tl *2 

= KF|2 X] E^^i ( S ) C i3,h ( U ) + (*2 - *l)Ci 2 ,i! (*2 - tl + M - S) 

+ CVjji^z - ti +u)Cj 3 ,t 1 (*2 - ti - s)} - C ll ^ 1 (s)C' l2 j 2 (u) + Q. 

So, \T n \cov{Ci ldl (s),d i2lh (u)} ^Y l rez{ G h,h( r ) C i2,ii( r + u - s ) + C i2,h {r+u)C hlil (r- 
s)} as n — * oo. Let T n (u) = {t : t e T„, t + u e T n }. We then have 

cov{C„ (hi ,Ui), C n (hj , ) } 
1 1 



|5(hj)| |S(h 7 - 



X! X! cov ] W~\ E Z{s k ,t k )Z{s k + \n,t k +Ui), 
ES(h i )s i es(h I ) U n ts,er„(Mj) 

1 *,eT„k) J 



|5(hi)| |S(h,-)| ^ ^ -i-*,* v-*y.-.,. v-j,j- 



Appendix B: Proofs of the theorems 



Proof of Theorem 1. Consider the covariance term, cov{C„(hi, m), C n (hj,Uj)}, where 
(hi,Ui), (hj,Uj) € A, and let D l n =D n (hi,Ui). Then, 

. ■ y^y^ COv{Z(s 1 ,ti)Z(Si +hi,ti +U i ),Z(s 2 ,t2)^(S2 +h 3 -,ti +Uj)} 

= ^ E E cov{Z(0,0)Z(h i , Mi ),Z(s,t)Z(s + h J , < + ,!,)} 

= E cov{z(o, o)z(h,, «o, z(s, t)z( s + h,-, t + Uj )} 1 J " n ^ ~ Jf; . 

Di-D' ' X ' 

Applying conditions (C2), (2.1) and Kroncckcr's lemma, we conclude that 

|Z?„|cov{C*„(h J ,u I ),5 n (h„u,)}^ E J2cov{Z(0,0)Z(h t ,u l ),Z(s,t)Z(s + h J ,t + u J )}. 

sgz 2 tez 
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Let a 2 = £ 8gZ2 £ tgZ cov{Z(0,0)Z(hi,ui),Z(s,t)Z(s + h 2 ,i + u 2 )} and A„ = 

y/\D n \{C n (h,u) - C n (h,u)}. By Lemma A.2, A n 4 N(0, cr 2 ). We apply the Cramer- 
Wold device to prove the joint normality. □ 

Proof of Theorem 2. See the proof in Li et al. (2007). □ 

Proof of Theorem 3. This is a special case of Lemma A. 4, setting p = 2 and q = 1. □ 

Proof of Theorem 4. Asymptotic normality of G„ follows immediately from Lemma 
A. 4, setting p = 3, q = 0. □ 
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